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PACS 67 . 85 . -d - Ultracold gases, trapped gases 

PACS 03.65.Vf - Phases: geometric; dynamic or topological 

PACS 05 . 45 . -a - Nonlinear dynamics and chaos 

Abstract. - We present relativistic linear stability equations (RLSE) for quasi-relativistic cold 
atoms in a honeycomb optical lattice. These equations are derived from first principles and provide 
a method for computing stabilities of arbitrary localized solutions of the nonlinear Dirac equation 
(NLDE), a relativistic generalization of the nonlinear Schrodinger equation. We present a variety 
of such localized solutions: skyrmions, solitons, vortices, and half-quantum vortices, and study 
their stabilities via the RLSE. When applied to a uniform background, our calculations reveal an 
experimentally observable effect in the form of Cherenkov radiation. Remarkably, the Berry phase 
from the bipartite structure of the honeycomb lattice induces a boson-fermion transmutation in 
the quasi-particle operator statistics. 
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Progress in condensed matter and particle physics has 
jbeen periodically marked by significant mutual exchanges 
between the two disciplines, many proposals for which 
are realized in model systems of ultracold quantum gases 
in optical lattices [T]-[3]. Recent active areas of research 
include holographic dualities such as AdS/CFT [4], the- 
oretical constructions of superstrings in ultracold quan- 
tum gases [5j , chiral confinement in quasi-relativistic Bose- 
Einstein condensates (BECs) and our own derivation 
of the nonlinear Dirac equation (NLDE) describing ultra- 
cold bosons in a honeycomb optical lattice f7| . Our inves- 
tigation into relativistic effects in BECs is motivated by 
this spirit of cross fertilization with the aim of tying in 
theory to experiment. 

In this Letter, we develop the relativistic linear stabil- 
ity equations (RLSE) for the NLDE. Moreover, we find 
emergent nonlinear localized solutions [8] to the NLDE, 
including solitons, vortices, skyrmions, and half-quantum 
vortices, the latter so- far unobserved in BECs. Although 
most of these objects have been studied in multicompo- 
nent BECs, such models lie within the usual Schrodinger 
many-body paradigm. In contrast to this paradigm, our in- 
vestigations reside within a relativistic framework in which 
the elementary excitations are governed by a Dirac-like 
equation. This provides a fundamentally different con- 
text distinguished by the presence of a Berry phase, so 



that exchange of two vortices leads to integer or half- 
integer exchange statistics. It is not surprising that the 
elementary excitations in our theory exhibit a rich struc- 
ture: a Dirac-like dispersion which obeys either bosonic 
or fermionic statistics depending on the strength of the 
contact interaction. Consequently, in order to determine 
the quasi-particle states and energies we cannot rely on 
the Bogoliubov-de Gennes equations (BdGE) since these 
are based on nonrelativistic quantum mechanics. Instead, 
we derive, from first principles, the RLSE which give the 
correct low energy dynamics for an arbitrary background 
condensate. The RLSE are reducible to the BdGE in cer- 
tain limits, and so may naturally be considered relativistic 
generalizations of the latter. Based on the RLSE we pre- 
dict Cherenkov radiation that can be measured in experi- 
ments: the combination of lattice and particle interactions 
results in a rich spatial distribution that is not seen in the 
BdGE for the uniform case [91. 

In the laboratory, the NLDE can be obtained by cool- 
ing bosons into the lowest Bloch band of a honeycomb 
optical lattice \iO\ ; the lattice is constructed by establish- 
ing three phase-locked interfering laser beams in a plane 
while freezing out excitations in the vertical direction as in 
Fig. [H To obtain the desired Dirac structure, particles are 
first condensed into the lowest energy state (zero crystal 
momentum) of the lattice and then adiabatically trans- 
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Fig. 1: The honeycomb optical lattice, (a) Cross section of 
the band structure showing K and K' points for gapped and 
ungapped systems, (b) The velocity and acceleration of the lat- 
tice, with A and B sub-lattices, are functions of the frequency 
offsets for interfering lasers, (c) Two-dimensional Dirac cone at 
K and K\ 

lated to the Dirac point at the band edge (see Fig. [T]) by 
adiabatically tuning the relative phases between the laser 
beams. We emphasize that the Dirac point, which is key to 
the NLDE and our predictions, is maintained even in the 
presence of the shallow harmonic trap endemic to atomic 
BECs [11] • Notably, the NLDE may also be obtained by 
the alternative method of using a square optical lattice 
with a staggered gauge field induced by a time-dependent 
optical potential ^ similar arrangement using only 
bosons is expected to show the same low energy struc- 
ture as in our model. Nonlinear phenomena in BECs have 
been studied extensively over the past decade |8J, from 
single-component vortices in rotating, trapped BECs [13] 
to complex mult i- component order parameters [l4l[T5] re- 
sulting from interactions between the different compo- 
nents and the possibility of nontrivial topological windings 
of the internal symmetry space around a singular vortex 
core. Some form of BdGE analysis plays a central role in 
such constructions as a means of probing stability as well 
as for gaining a deeper understanding of the low-energy 
fluctuations. 

Derivation of the RLSE. — Since the system we 
describe in this letter is a BEC confined strictly to two 
spatial dimensions, it is appropriate to recall the justifica- 
tion for such a construction before presenting the RLSE. 
For uniform 2D systems the Mermin- Wagner theorem for- 
bids the formation of a true condensate defined by an in- 
finite phase coherence length. This comes from the fact 
that the density of states diverges in the 2D case for finite 
T. Instead, one sees the formation of a quasi- condensate 
characterized by local phase coherence restricted to finite 
size regions. The size of these regions greatly exceeds the 
healing length so that all of our solutions are realizable 
in this picture. However, the inclusion of a harmonic con- 
fining potential allows the formation of a true 2D conden- 
sate. The potential places a lower bound on the energy for 
fluctuations and, since it is these long wavelength fluctua- 



tions that are responsible for destroying long range order, 

the trap provides a means of expanding the spatial range 

of vahdity of the mean field description. 

To obtain the low energy excitations of solutions of 

the NLDE, we must find the correct set of equations 
that describe quasi-particle states analogous to the BdGE 
equations for the general case. These are obtained from 
the Hamiltonian for a weakly interacting Bose gas, H = 
JdriP'^Ho^P^^ Jdr^P^^iP, Ho = h^\/^ /2M ^V{r), and 
working through four steps |16J. (1) Take = ^c(r) + 
J?/^q(r) (condensate + quasi-particles), with 61^^ small. (2) 
Impose a constraint on to eliminate linear terms in 
(5?/;q(r), keep only quadratic terms in (5V;q(r), and expand 
as a sum of particle and hole creation operators. (3) In- 
voke Bloch-state expansions for ^c(r) and (5?/^q(r) and take 
the lowest band. (4) Take the long-wavelength limit while 
taking momentum with respect to the Dirac point K, such 
that k <C q <C K, where q is the momentum of the con- 
densate relative to the Dirac point K and k is the quasi- 
particle momentum measured relative to q, and finally di- 
agonalize the quasi-particle part of the Hamiltonian. One 

finds the RLSE: . 

^Uk - /7^Vk = ^kUk, (1) 

#*Vk - t/^Uk = -^kVk , (2) 



where the matrix coefficients are defined as 

§ = diag(|^A|',|^B|'), (3) 
^k = diag(^k,^k), (4) 
= meff-/i + 2/7|^A|'-^V(/)A- V 

+ |V(/)A|-i(V2(/>A), (5) 

[#]2^2 = meff - /i + 2/7 \^b\^ - iV(t)B ■ V 

+ |V(/>b|-^(V2(/>b), (6) 
V\ (7) 
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J2,l 



Here, V = {dx + idy) is the single particle Dirac op- 
erator. Also, ^ = is the BEC order param- 
eter at the K Dirac point, with normalization on sub- 
lattice components J drd^^P + I^Bp) = 1- Analogous 
equations hold for the inequivalent Dirac point at — K. 
Cast in this highly compact form, Eqs. ([I])-© are rem- 
iniscent of the BdGE and may be solved for the spinor 
quasi-particle amplitudes Uk(r) = [^ik,A(r), iik,B(r)]"^ and 
Vk(r) = ['^k,A(r), '^k,B(r)]"^ and the quasi-particle energy 
Ek. The components of these 2-spinors represent quan- 
tum fluctuations of the sublattice condensate order pa- 
rameters and ^ B which in general are nonuniform C- 
functions on the plane. The presence of the local phase 
of the condensate (/>A(B)(r) indicates the complex inter- 
action between the local superfluid velocity of the con- 
densate V5^^(^)(r) = V0^(B)(r) and the spinor quasi- 
particles Uk(r) and Vk(r). We have taken h = ci = 1 
for simplicity, where ci is the effective speed of light in the 
NLDE. Note also that we have included an effective mass 
meff (anisotropic lattice) that competes with the chemical 
potential /i. 
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The RLSE exhibit a positive-negative energy inversion 
symmetry which is found by complex conjugating Eqs. ([T]) 
and (j2j) followed by a spatial parity inversion. This is simi- 
lar to the BdGE where negative energy modes are allowed 
by symmetry but are ignored since the underlying Hamil- 
tonian is positive-definite. It is important to note that for 
a moving condensate, the negative-energy modes cannot 
be removed and are crucial indicators of Cherenkov ra- 
diation. However, in our case, the Dirac Hamiltonian is 
not positive-definite since our theory is defined at zero 
lattice energy, not the lowest energy Bloch state, so we 
must respect the presence of both energy raising and low- 
ering modes. Another important feature is that the RLSE 
are reducible to the BdGE when the local lattice potential 
energy is the main contributor to the condensate chemi- 
cal potential and the condensate is slowly varying (quasi- 
uniform background), i.e., >> /7, Ej ^ |Eo|, 
where Eq = T^oa{b) - jdYw\^-^^HQWA{B) is the lo- 
cal self energy for an arbitrary lattice site, with Wa{b) = 
w{t — r^(B)) the Wannier functions. 

Quasi-particle Operator Statistics. — Since we 
use only bosons in our construction, and we have shown 
that these collectively produce the spinor structure of 
Dirac theory, it is natural to ask what are the quantum 
statistics for collective excitations. When quantum fiuc- 
tuations are weak, quasi-particle states are superpositions 
of direct products of bosonic operators and spinor single- 
particle states: operators retain their bosonic structure. 
In contrast, when the depletion of the condensate becomes 
significant, quasi-particle operators at each lattice site be- 
come sums of superpositions of terms each with a different 
phase, i.e., the single-particle bosonic operators b and 
become entangled with the on-site phase. We can illus- 
trate this more rigorously. Consider a BEG in the thermo- 
dynamic limit: and N/A4 = constant, where N 
is the total number of particles and Ai is the number of 
lattice sites. Also, we focus on the regime where U/th <C 1 
at T = 0. The ground state is then described by the su- 
perfiuid density order parameter and the quasi-particles 
are coherent long wavelength fiuct nations in the phase. A 
Gutzwiller ansatz in terms of on-site number states pro- 
vides an adequate formulation of the wavefunction. If we 
consider increasing the particle interaction U towards the 
critical value Uc that separates the superfiuid from the 
Mott-insulating phase, we would expect to see an increase 
in depletion of the condensate. The part of the wavefunc- 
tion that describes particles outside the condensate has 
complete phase decoherence so that the phase fiuctuations 
at the i^^ site, for N- particles outside the condensate, can 
be described using the eigenstates of the phase operator 
in the number state basis [T7]: 
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where |0)- is the ground state wavefunction for the i^^ 
site. With this definition, the i^^ site phase operator is 



is a reference angle which must be averaged over when 
computing observables. In general, the depletion number 
Nl is a function of the total number of particles Ni and 
the interaction U. When N- is large, this dependence on 
N and U can be incorporated into a coefficient fi{N^U) 
that multiplies each term of the sum while maintaining 
the relative sizes Ni << N- « 1: 



n=0 



|0), = Dia) |0), . 
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where = /ie*^° and we have used the symmetric form for 
the phase coherent state D{a) = ^'^i^.-a-bi^ r^^^ normal- 
ized quasi-particle operators can now be defined having 
the correct Bogoliubov limit: 

al = bl D\a)/^/lTW, a, = D{a) k/V^ + H^ ■ (10) 
With this definition, the average number of bosons form- 
ing a quasi-particle is (h^b^ = l + where \a\ en- 
codes the degree to which the bare quasi-particle oper- 
ators b and b^ are dressed or fused to the background 
phase. When working at the lattice scale it is clear that 
the reference angle 60 must be single-valued under one 
full rotation so that the dependence on the polar angle 
is e*^^°+^^. However, when translating to the continuum 
limit, the two triangular sublattices of the honeycomb lat- 
tice become identified with the same spatial point so that 
a 2 — 1 mapping of parameters is needed consistent with 
the double covering map p : SU{2) ^ S0{3). The contin- 
uum quasi-particle statistics can then be determined by 
computing the Berry phase (holonomy) by adiabatically 
transporting a quasi-particle through a suitable closed 
path and correctly accounting for the SU(2) structure of 
the single-particle states. The path must remain within a 
degenerate subspace of the Hamiltonian so that we may 
isolate the geometric phase from the dynamical phase. 
The holonomy is 7 = exp §^ dO {^a{0)\ I'^aiO)) with 



|0) and a = /e^(^o+^)/2. A straight 



forward calculation yields 7 = exp(z7r|ap) demonstrating 
fermion exchange statistics when |ap = 1. 

Physical Parameters and Regimes. — We list 
first the fundamental dimensionful parameters that we 
use. They are as follows: the average particle density 
no, the chemical potential /i, the lattice spacing a, the 
s-wave scattering length a^, the mass of the constituent 
bosons M, and the lattice well depth Vq. Several rele- 
vant composite quantities may be constructed from these. 
These are the effective speed of light ci = thCiV^/2h^ 
the sound speed Cs = ^s/Uno/M, the interaction strength 
U = Anhas/M^ the healing length ^ = thCi\/?t / 2hn{)U ^ and 
the hopping energy th = jd^r w"" Hqw^ where th depends 
on a and Vq^ respectively, through the overlap of Wannier 
functions and the lattice potential inside Hq. Two fun- 
damentally important constraints regarding these quan- 
tities should be stated. First, in order to avoid reaching 
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the Landau velocity at the band edge and creating un- 
wanted excitations we require ci < c^, where Cs is the 
sound speed. Thus we require thay/3/2h < y^Uno/M or 
{thaV3/2h){47rh^asn/M)-^/'^ < 1. For ^'^Rb with th = 
h X 10^ Hz, a = 0.5 X 10"'^ m, as = 5 x 10~^m, n = 
2 X lO^^cm"^, we get ci/cs ^ 0.17. Second, in order for 
our long- wavelength approximation to be correct, we re- 
quire the NLDE healing length ^ = thaV3/2hnoU a; 
using the same values for the physical parameters, we find 
^ ^ 9.43 a. 

Next, we discuss the physical regimes for our theory. 
First, consider the extreme weakly interacting regime, i.e., 
noU/th <C 1. Then excitations of the ground state obey 
bosonic statistics. Furthermore, at length scales much 
larger than the healing length, ^/c <C 1 where /c is a char- 
acteristic quasi-particle momentum, excitations are com- 
prised of correlated particle-hole pairs that propagate with 
a dispersion given by £^ oc /c^/^. This is a Bose gas of com- 
posite particles in the sense that excitations of opposite 
spin are paired up (albeit non-locally) to form bosons. In 
contrast, for ^/c ^ 1, excitations are particle-like which 
corresponds to the case where spin eigenstates are excited 
independently. These states reflect the bipartite structure 
of the lattice (multi-component) but are local objects and 
so also reflect the bosonic nature of the fundamental con- 
stituent particles. They exhibit a mixture of fermionic and 
bosonic properties having a Dirac-like dispersion ex k but 
with quantum operators that obey bosonic statistics. This 
is a hybrid Dirac-Bose gas. On the other hand for mod- 
erate interactions, noU/th 7^ and close to the critical 
point on the superfluid side, higher order processes be- 
come significant and quasi-particles become highly non- 
local objects. Here quasi-particles are heavily "dressed" 
in the presence of the background condensate and so ac- 
quire a geometric phase resulting in anticommutation re- 
lations for their associated creation and annihilation op- 
erators (see previous section). In this regime, the Bogoli- 
ubov approximation fails to provide a good model but we 
expect that the Dirac spinor structure will remain robust 
against many-body effects with quasi-particles undergoing 
a renormalization due to many-body interactions similar 
to the case of graphene [18] . This describes an interacting 
Dirac gas. 

Uniformly Moving Condensate. — Now we return 
to the RLSE and solve them for the simplest case of a uni- 
form background ^(r) = y^e*^"'^(l, Cq)^, where Co G C 
contains a relative phase, no is the average particle den- 
sity, and q is the condensate momentum measured with 
respect to the Dirac point. In order to obtain the coher- 
ence factors and quasi-particle dispersion, we must then 
solve a 4 X 4 eigenvalue problem; the RLSE yield 

£;k = cihq ■ k ± ^{cihkY + noUcihk . (11) 

In keeping with the usual Bogoliubov notation found in 
the literature, we may write = (cJ/Q)q ■ ^ 
where = cifik is the single quasi-particle energy for 



zero interaction and = \/ (c^)^ + ^oU is the quasi- 
particle energy for a static background. The associated 
coherence factors can then be written as \uj^^a{b)\ = 

{El + Cihk)/^AElcihk, \Vk^A{B)\ = |^fc,A(B)|(+ ^ -). 

The full interacting Hamiltonian is given by ^rlse = 
\UnlA + cihq - Y.'k(^4 + ^o/7) + Y.'k ^^4^^ ' where A 
is the area of the plane. The first three terms are the 
mean-field and quantum corrections to the condensate 
energy and the last accounts for the number of quasi- 
particles present in the system. The constant c[ is de- 
fined in terms of the overlap integral between Wannier 
states at neighboring lattice sites by c5 = \f?>aT/2h^ where 
TA,B - Jdrw\VwB and r = |rA,B|u. 

The low energy behavior of the uniform condensate has 
a rich structure. The q = case corresponds to a con- 
densate with zero crystal momentum measured from the 
Dirac point but with momentum K relative to the low- 
est Bloch state of the crystal. The idea of a condensate 
in motion relative to its background has been treated in 
both free-space as well as the case of a moving background 
lattice [T9l[20]. Physically, the lattice potential is moving 
relative to the stationary condensate (laboratory frame). 
Two-body collisions reduce the momentum of some parti- 
cles relative to the lattice (slowing down) and increase the 
momentum of others (speeding up) corresponding to a fi- 
nite depletion of the condensate. In the laboratory frame, 
a two-particle collision appears as one particle gaining a 
component of momentum to the left and the other a com- 
ponent to the right. This is consistent with the well known 
particle-hole symmetry of the Dirac Hamiltonian: neg- 
ative energy states can be interpreted as positive energy 
states that propagate in the opposite direction. In our the- 
ory these are quasi-particles with momentum K — k (for 
the K-Dirac point) relative to the lowest Bloch state. 

For q = then, we get E^"^^ ±E'l = 

(^fe)^ + noUe^. The two energy regimes evident here 
are separated by the condition cihk/noU = £,k ^ 1. At 
short wavelength, k(_ ^ 1 so that E^^ ^ ±{cihk-\-noU /2) ^ 
where the dominant first term reflects only the presence 
of the honeycomb lattice, while the second term is a small 
mean-field Hartree shift due to the interaction with the 
background. When /c^ <C 1, we find E^^^ ^ ±y^/c/^. 
These are collective excitations induced by the particle 
interactions just above the condensate energy. The pres- 
ence of negative energy modes means that the conden- 
sate may lower its energy through spontaneous emission 
of radiation. This process can be suppressed by intro- 
ducing an anisotropy in the lattice by breaking the A- 
B sublattice degeneracy with a deeper optical lattice in 
one direction [2P. This results in an additional term 
in the dispersion opening up a mass gap 2meff at the 
Dirac point. For the negative energy modes we then 
have Ejr\mef^) = 2meff — \/ (e^)^ + rioUe^ so that ex- 

■"^Note that whereas for the effective speed of Hght we have [ci] = 
m- , in contrast [cj] = ■ since ta^b is an integral over the 
gradient operator rather than the Laplacian. 
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citations require a minimum momentum determined by 
cihkmin = \/4m^^ + nQ^/^. Alternatively, we can con- 
sider the effect of the confining potential: this sets a 
lower bound for quasi-particle energy given by l^^^^j^i^)! ^ 
^J (cih2i\ I -\-noU Cih27r/R_\_, where R± is the char- 
acteristic trap radius in the 2D plane. 

Cherenkov Radiation. — This usually refers to the 
anisotropic emission of electromagnetic radiation from a 
source whose speed exceeds the local speed of light in some 
medium [22] . This concept generalizes to any source mov- 
ing through a medium at a speed that exceeds the phase 
velocity of the elementary excitations of the medium. For 
example, a BEC moving in the laboratory frame, or with 
respect to a background, will "radiate" (emit particles) 
when its speed exceeds the sound speed. Moreover, the 
radiation will be emitted in a cone subtended by a specific 
angle in the direction opposite the motion of the BEC. 
The RLSE can be used to demonstrate this effect in the 
present context of a BEC in a honeycomb optical lattice. 

For a BEC with momentum q > measured from the 
Dirac point, examination of the angular dependence of 
Ek reveals an intriguing structure for the emission of 
Cherenkov radiation. We observe the following properties 
for Ek. (1) When v < ci^ where v = c^q is the condensate 
speed, all excitations have positive energy regardless of the 
angle of emission. (2) When v > ci^ quasi-particle energies 
are positive only for emission angles (measured relative to 
q) for which < Oc = cos~^ {—ci/v) while all other modes 
have negative energy corresponding to the emission of ra- 
diation in a backwards cone bounded by Oc. When v = ci^ 
Oc = 7T marks the onset of radiation, in which case ra- 
diation is only emitted in the direction opposite q. This 
unique directional property of the radiation suggests an 
obvious detectable signature in the laboratory: a time-of- 
fiight analysis of a BEC prepared with precise values of 
the parameters should show a predictable shift in the mo- 
mentum distributions between the forward and backward 
directions. 

Nonlinear Localized Modes. — To obtain solutions 
of the NLDE that are localized in x, y for U> 0, we substi- 
tute the plane-polar ansatz ^^(r) = exp[zp^(^)]FA(r), 
^^(r) = CBexp[ipB{0)]FB{r) into the NLDE. Then ca = 
i, cb = 1, and there are two possible combinations for 
the angular functions: (i) pa{0) = {I — 1)0^pb{0) = 
10 ; (ii) pa{0) = {I - 1/2)0, pb{0) = (/ + l/2)i9, with 
/ G Z. In particular, / = in (i) corresponds to a 
vortex configuration in ^ a filled in at the core with 
a nonzero soliton for ^b- Solutions of this type exist 
for different relative values of /i and U and for several 
asymptotic values of the components: lim^^oo(^A, ^ b) ^ 
{(-iV^,0),(-iv^, V^),(0,0)}. For / = 1, we obtain 
the same types of solutions but with ^a and ^b ex- 
changed. For / > 1, centripetal terms are present for both 
FA{r) and ^^(r) so that we must have ^a(O) = ^b(O) = 
and both components are vortices with zero core den- 




Fig. 2: Localized solutions of the NLDE, (red dashed 
curves) and (blue solid curves), in units of NLDE heal- 
ing length ^: (a) vortex/soliton, (b) ring-vortex/soliton, (c) 
half-quantum vortex, (d) planar skyrmion, (e) line skyrmion, 
and (f) line-soliton. 

sities. For the / = 1 case, we also obtain a skyrmion 
solution for which the pseudospin S = ^(r)cr^(r) (with 
Pauli vector a) exhibits an integral number of fiips near 
the core and approaches a constant value far from the core. 
This feature is encoded in a topologically conserved charge 
(I/Stt) J^dre^^S • diS x djS which one recognizes as the 
Pontryagin index that classifies the mapping S^^^^ Sqq 
where the two circles S^^^^ and Sqq parameterize the ro- 
tations between the densities Pa{b) and the polar angle 
rotation on the 2D boundary dQ at spatial infinity. In 
general, similar types of solutions exist for (ii) above. An- 
alytical and numerical solutions are plotted in Fig. [2] for 
which different values oi ji/U and / allow us to obtain the 
different asymptotic forms. 

Besides vortices with integer phase winding, we also 
find solutions with fractional phase winding, called half- 
quantum vortices (HQVs). Ordinarily, analyticity (single- 
valuedness) of the order parameter forbids the rotation 
of the phase of ^ around the core to take on fractional 
values. In the NLDE, ^ can acquire a coherent internal 
Berry phase in addition to an external phase whose an- 
gles are identified with the polar angle [23l[24]. Such 
states may have half-integer winding in both the inter- 
nal and external phase angles while remaining single- 
valued overall. We obtain HQVs with asymptotic form 
lim^^oo ^HQv(r) = 2zy^e-^^/2[cos(^/2),zsin(^/2)]T; 
the complete solution is shown in Fig. [2] (c). 

We also obtain one-dimensional kink-soliton, skyrmion, 
and line-soliton solutions. The kink and skyrmion solu- 
tions are obtained by a straight-forward substitution of 
the ansatz ^{x) = r] [cos{(p) , s'm{(p)]^ into the NLDE and 
then considering the distinct cases where (p = constant 
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Fig. 3: Plots of lowest quasi-particle excitation for the vor- 
tex/ soliton configuration. 

(kink) 017] = constant (skyrmion). The line-soliton so- 
lution is obtained when both 77 and if are functions of x 
with the additional condition that, at the origin, r] remains 
below a certain value. This ensures that fi^ < \/U/^ < 
{U + 1)/2U ^ U < 3.365 and /i^ < 0.649, which allows 
the wavefunction to collapse away from the 7/-axis while 
the nonzero wavefunction near and along the y-ax.is has 
a Lorentzian form in the x-direction due to the attractive 
effect from the kinetic terms. 

Localized Mode Stability. — To check the stability 
of our localized solutions we substitute the correspond- 
ing solution into the RLSE and obtain the low energy 
spectrum. For the case of the vortex/soliton (/ = 0), 
we solve the eigenvalue problem numerically and obtain 
Eo = -3.9274 + 0.0020 i in units of noU. Similar to the case 
of a trapped rotating BEC [25], we find that the lowest 
mode in the spectrum is anomalous with negative energy 
and positive norm. The coherence factors shown in Fig. [3l 
near the core of the vortex, appear as |i^/c,a(b)P ^ 10~^, 
b/c,A(B)P~10"^ SO that Uk^AiB) » Vk,AiB)' This mode is 
dynamically unstable due to the presence of a small imag- 
inary component which gives a decay rate relative to char- 
acteristic oscillation time of 0.0020/3.9274 = 0.0005. The 
lowest quasi-particle energies for the other localized solu- 
tions are: -3.9276 + 0.0019 i; 2.634 x 10^ + 9.96 x lOH; 
-3.9274 + 0.0019 i; 7.8409 x 10"^- 9.9993 x 10^ z; 7.9349 x 
10"^- 9.9993 X 10^ i; for the ring- vortex/soliton, half- 
quantum vortex, planar skyrmion, line-skyrmion, and line- 
soliton respectively, where all quantities are given in units 
of rioU. 

In conclusion, we have shown that an effective relativis- 
tic fermionic system may be designed using ordinary cold 
bosonic atoms as the underlying degrees of freedom. We 
solved the resulting NLDE for different classes of nonlin- 
ear modes including half-quantum vortices. We derived 
and solved relativistic linear stability equations and gave 
explicit criteria for experimental observation of Cherenkov 
radiation, as well as predicting an anomalous mode for the 
vortex/soliton solution. Density profiles may be observed 
by time-of-flight techniques to detect both massive and 



massless Dirac fermions in the laboratory [2T1[26]; non- 
linear modes involving phase winding can be created by 
techniques analogous to those used at JILA [2T; and we 
anticipate that Bragg scattering can be used to populate 
the Dirac cones at both K and points, leading to ar- 
bitrary superpositions over our localized solution types 
between the two cones, and thereby populating all four 
components of the Dirac spinor. 
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